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field source detection: ROSAT sources in NGC 6397 
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ABSTRACT 

We present a two-dimensional version of the classical one-dimension al Kolmogorov- 
Smirnov (K-S) test, ex tending an earlier idea due to Peacock (1983) and an imple- 
mentation proposed by Fasano fc Franceschini (1987 ). The two-dimensional K-S test 
is used to optimise the goodness of fit in an iterative source-detection scheme for 
astronomical images. The method is applied to a ROSAT x-ray image of the 
post core-collapse globular cluster NGC 6397 to determine the most probable source 
distribution in the cluster core. Comparisons to other widely-used source detection 
methods, and to a Chandra image of the same field, show that our iteration scheme is 
superior in measuring statistics-limited sources in severely crowded fields. 

Key words: methods: data analysis - methods: statistical - globular clusters: indi- 
vidual (NGC 6397) - x-ray: stars 



1 INTRODUCTION 

Deep x-ray imaging of crowded fields, even with increas- 
ing angular resolution and sensitivity (with Einstein and 
ROSAT, and now Chandra), is invariably limited by the 
small number of source counts and by the relative size of 
the point-spread function (PSF) compared to the angular 
separation between the objects. Determining the underlying 
source muflguitiLiuu In sudi a itiglmy Is often byyoud thy 

capabil itica of conventional aourec - finding algoritlinia. 

Classical x-ray source detection methods are based 
on a sliding detection cell of a fixed size across the im- 
age, and calculating the signal-to-noise ratio (S/N) at each 
step. To find S/N, common detection algorithms for pro- 
cessing data from Emstetn and ROSAT (implemented in 
IRAF Q/PROsQ) use either an average background (as deter- 
mined from a source-free section of the image) or a local 
background (from a region around the detection cell). How- 
ever, both methods fail to discern blended faint sources in 
crowded fields where the background is affected by over- 
lapping PSFs. Source detection is somewhat improved by 
image deconvolution, e.g. with the Lucy-Richardson (L-R) 
algorithm or with the Maximum Entropy Method (MEM), 
or by wavelet smoothing. Deconvolution algorithms pro- 
vide higher positional sensitivity in moderately crowded 



^ Image Reduction and Analysis Facility; developed and main- 
tained by the National Optical Astronomy Observatories. 
^ Post- Reduction Off-line Software; developed and maintained by 
the Smithsonian Astrophysical Observatory. 



fields, but suffer from such undesirable effects as noise- 
amplification and "leakage" (associating counts from fainter 
sources to brighter nearby ones). Wavelet detection imple- 
mented as task wav detect in the Chandra pro cessing pack- 



age (available at http://asc.harvard.edu/ciac) does well in 
crowded fields, provided the sources are either sufficiently 
separated (>3-5 FWHM, or >3-5 arcsec) or within ~2- 



3 FWHM and a re similar in flux ( Damiani et al. 1997 ; Free- 



man et al. 2001). When the PSFs are heavily blended (sepa- 



ration between the source centroids < 1.5 FWHM), individ- 
ual sources cannot be distinguished and their relative fluxes 
cannot be measured. A superior source-detection method is 
needed for severely crowded fields containing multiple fain t 
sources, e.g. globular cluster cores (Hertz & Grindlay 1983), 
or nuclear bulges in external galaxies. 

A powerful technique to compare statistics- limited sam- 
ples is the Kolmogorov-Smirnov (K-S) test which, unlike its 
alternative - the Pearson test - does not require binning 
of the data. Unfortunately, the classical K-S test is applica- 
ble only to one-dimensional distributions, and any attempts 
to convert a two-dimensional image to one dimension (e.g. by 
collapsing it onto a vector, or by azimuthal binning around 
a point) lead to unwanted loss of information and power. 
For some time now, a multi-dimensional version of the K- 



S test has been known (Peacock 1983; Gosset 1987), which 
performs better than the y test in the small-number statis- 



tics case (Fasano fc Franceschini 1987, hereafter, FF), and 



can be successfully applied in parameter point estimation 
in a manner similar to the widely used maximum-likelihood 
(ML) method. These properties of the multi-dimensional K- 
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S test make it viable for incorporation in source-detection 
methods. 

In this paper, we re-visit the characteristics of the K-S 
test in two (and three) dimensions and examine its power in 
comparing different realisations of crowded low S/N fields. 
As an application of the test, we devise an iterative source- 
modelling scheme that aims to minimise the K-S statistic in 
search of the optimum underlying source distribution in an 
image. Based on our Monte Carlo simulations, we find that 
our iterative algorithm is a powerful tool for faint object 
searches in crowded fields. We apply the algorithm to deter- 
mine the faint x-ray source distribution in a deep ROSAT ex- 
posure of the post core-collapse globular cluster NGC 6397 
which has also been an alysed with ML techniques by 



bunt I 

tions with those of fVerbunt & Johnston 



Ver- 



Johnston (200C ). We compare th e derived x-ray P' 



Grindlay et al. 2001a) identifications 



quent optical (HST, Taylor et al. 2001) and x-ray {Chan 



and with our subse- 



2 THE TWO-DIMENSIONAL K-S TEST 
2.1 Description 

The classical one-dimensional (ID) K-S test makes use of the 
probability distribution of the quantity Dks, defined as the 
largest absolute difference between the cumulative frequency 
distributions of the parent population and that of an n- 
point sample extracted from it. Since Dks is approximately 
proportional to l/yTi, one usually refers to the probability 
distribution of the quantity Z„ = DKSy/n. For a given n, 
the values of Zn corresponding to a given significance level 
SL (denoted as Z„^sl) increase slightly with n. For large 
n, the integral probability distribution P(> Zn) = 1 — S L 



appr oaches the asymptotic expression (Kendall & Stuart 



1970 



P(> Zn) = 2 ^(-1)*-^ exp(-2fc2z„) 
fc=i 

which is satisfactory for n ^ 80. For the two-sample K-S test, 
which compares distributions of different sizes (ni and n2), 
the probability distribution P{> Z„) remains unchanged 



provided that n is set to 



The ID nature of the test 



implies that it does not depend in any way on the shape of 
the parent distribution. 

In a two-dimensional (2D) distribution, each data point 
is characterised by a pair of values, {x,y). As with the ID 
K-S test, the maximum cumulative difference between two 
2D distributions is found over the {x, 3/)-plane. In the case 
of distributions in more than one dimension however, the 
procedure to cumulate the information onto the plane is not 
unique. FF made use of the total number of points in each of 
the four quadrants around a given point {xi,yi), namely, the 
fraction of data points in the regions {x < Xi,y < yi), (x < 
Xt,y > yi),{x > Xi,y < yi),ix > Xi,y > yi). The 2D statis- 
tic Dks is defined as the maximum difference between data 
fractions in any two matching quadrants of the sample and of 
the parent population, ranging over all data points. The Zn 
statistic is defined similarly as in the ID case, Zn = DksV^, 
where for a two-sample 2D K-S test n = "'-"^ . 

Based on their Monte Carlo simulations, FF deduce 
that the 2D integral probability distribution P(> Zn) de- 



pends solely on the correlation coefficient (C'C) of the model 
distribution, i.e. that for a given CC the distribution of Zn 
in the 2D K-S test is (nearly) independent of the shape of the 
model, as in the classical ID K-S test. FF also observe that 
in the two-sample case, it is sufficient to take the average of 
the correlation coefficients CCi and CC2 of the samples as 

an estimate of CC. 

An important distinction between Peacock's and FF's 
version of the 2D K-S test was pointed out to us by the 
referee, which makes the latter generally less stringent. FF 
restrict the search for the maximum cumulative difference 
Dks to loci harbouring a data point, thereby often missing 
the location of the true maximum difference, which is almost 
always found for {x < Xi,y < yj), where i and j are two dif- 
ferent data points. Nevertheless, the maximum cumulative 
difference computed in such a way will have a tendency to 
vary in the same manner as the true maximum difference. 
Thus the FF statistic is probably well-behaved, at least as 
long as the genuine parent popula tion distribut ion and the 
assumed one are not too different (Cosset 1987). The latter 
situation is indeed expected when comparing distributions 
of point-spread functions in two images. The advantage of 
FF's approach is speed: order n, instead of n^. The disadvan- 
tage is its approximate nature and that the Dks statistic 
is sensitive to the correlation coefficient CC of the distri- 
butions, requiring its inclusion as a free parameter in the 
reference tables. Our Monte Carlo experiments below take 
into account both factors. 



2.2 Monte Carlo experiments 

Following the procedure in FF, we used the 2D K-S test corn - 
puter code provided in Numerical Recipes ( Press et al. 1997 ) 
to run our own Monte Carlo experiments. We studied the Z„ 
statistic by means of a Monte Carlo procedure using a uni- 
form distribution (CC = 0) within the unit square as a par- 
ent population. The analysis comprised of cases with number 
of points n per sample ranging from n = 5 to n = 50000. 
For any given n we produced a large number of simulations 
(from 100000 for n = 5 to 1000 for n = 50000), enabling us 
to construct the integral probability distribution P{> Zn) 
with sufficient accuracy. Values oi Zn,SL for uniform samples 
of all tested sample sizes n are listed in Table |l| 

Figure ^ presents a comparison between our and FF's 
results for the critical values Zn,SL as a function of n. Both 
sets of Monte Carlo simulations show similar tendencies in 
the behaviour of the Zn statistic, and are indistinguishable 
from each other for n ^ 500 within the statistical uncertain- 
ties (« 5 per cent for high n, due to the limited number of 
simulations). There is however a marked inconsistency be- 
tween the two sets of data in the low-SL, small-n part of the 
graph, where at the 30 per cent significance level, the values 
of Zn,SL differ by a factor of ~ 1.15. The difference is highly 
significant, given the fact that both FF's and our results 
for small n are based on 100000 simulations. We attribute 
this discrepancy to a detail in the implementation of the 2D 
K-S test: in particular, to whether the data point around 
which the Dks statistic is computed, is included in one of 
the cumulative quadrants or not. This discrepancy becomes 
irrelevant, when the same implementation of the 2D K-S test 
is used in both building the reference tables for the critical 
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Table 1. Critical values Z„ gi^ for a uniform uncorrelated distribution. 





SL(%)'^ 


30 


40 


50 


60 


70 


80 


90 


95 


99 




# of simul. 




















5 


100000 


0.76 


0.82 


0.88 


0.95 


1.03 


1.13 


1.27 


1.38 


1.62 


10 


100000 


0.79 


0.85 


0.92 


0.98 


1.06 


1.15 


1.29 


1.41 


1.65 


20 


100000 


0.85 


0.91 


0.97 


1.03 


1.10 


1.20 


1.33 


1.45 


1.69 


50 


100000 


0.91 


0.97 


1.04 


1.10 


1.18 


1.27 


1.41 


1.53 


1.78 


100 


10000 


0.96 


1.02 


1.08 


1.15 


1.22 


1.32 


1.45 


1.57 


1.81 


200 


5000 


1.01 


1.07 


1.13 


1.19 


1.27 


1.37 


1.51 


1.62 


1.86 


500 


5000 


1.06 


1.12 


1.18 


1.25 


1.33 


1.43 


1.57 


1.68 


1.91 


1000 


5000 


1.09 


1.15 


1.21 


1.28 


1.35 


1.45 


1.60 


1.73 


1.98 


2000 


5000 


1.11 


1.17 


1.24 


1.30 


1.38 


1.47 


1.61 


1.74 


1.99 


5000 


2000 


1.14 


1.21 


1.27 


1.34 


1.42 


1.50 


1.64 


1.74 


1.96 


10000 


1000 


1.15 


1.22 


1.28 


1.35 


1.43 


1.54 


1.69 


1.82 


2.13 


20000 


1000 


1.14 


1.21 


1.29 


1.37 


1.45 


1.56 


1.74 


1.89 


2.20 


50000 


1000 


1.16 


1.24 


1.32 


1.40 


1.50 


1.60 


1.78 


1.91 


2.19 


^Significance level {SL 


= 1 - 


P(> ^n)). 















tSize of sample. 




Figure 1. Critical values of the statistic Z„ as a function of 
sample size n, for values of the significance level SL varying from 
30 per cent to 99 per cent. The continuous lines are based on data 
listed in Table ^ and the dashed lines are based on data reported 
in Table Al of FF. 



values of Z„, and in comparing actual distributions using 
these tables. 



3 USING THE TWO-DIMENSIONAL K-S 
TEST ON ASTRONOMICAL IMAGES 

FF applied a 2D K-S test to astronomical distributions (not 
images) and showed that it can reject wrong hypotheses at 
a much higher significance level than the test. We extend 
the application of the 2D K-S test to comparing images of 
crowded fields, where it is of great interest to determine 
whether a proposed source distribution corresponds to the 
observed one. When the PSFs of the individual sources are 
heavily blended (source separation < 1.5 FWHM), classical 
x-ray source-detection methods fail to distinguish the indi- 
vidual obj ects. Crowded - field o ptical photometry tools (e.g. 
DAOPHOT, Stetson 1987, 1991) are also not suitable for the 



small-number Poisson statistics of x-ray images. Even the 
recently int roduced detection al g orithm based on wa velet 
transforms (Freeman et al. 2001; Damiani et al. 1997, and 



analysis package (CIAO^, does not produce adequate re- 
sults in the regime of severe source confusion and small num- 
ber of counts per source. 

Below we describe an implementation of the 2D K-S test 
to astronomical images. Provided that the PSF of the (un- 
resolved) sources in a crowded field is known and constant 
over the image, by comparing the image to a simulation of a 
proposed source distribution, we can obtain the K-S proba- 
bility P(> Zn) that the image and the simulation represent 
the same source configuration. The obtained probability can 
be used as a measure of the accuracy of both the positions 
and the intensities of the proposed sources, as well as an in- 
dication of the necessity for additional sources to account for 
the photon distribution. Because of the high sensitivity of 
the 2D K-S test, we expect that it should be able to discern 
positional discrepancies equal to a fraction of the FWHM of 
the PSF. 



3.1 Implementation of the two-dimensional K-S 

test 

3.1.1 A two-dimensional vs. a three-dimensional K-S test 

Astronomical images (e.g. from CCDs) have three dimen- 
sions: X and y pixel coordinates, and pixel intensity. Al- 
though the spatial distribution of photons that strike the 
detector is two-dimensional (two photons never fall at the 
exact same position) , they are binned by the detector in in- 
teger bins corresponding to the digital pixel size. Thus the 
incoming 2D distribution of photons with real- valued coordi- 
nates is transformed into a three-dimensional (3D) one with 
integer coordinates. It is possible to apply the 2D K-S test 
to an integer- valued 3D distribution by simply re-calculating 
the two-dimensional Dks statistic for every count in a given 
pixel. However, since the K-S test is designed to deal with 
real-valued distributions, any point in the {x, y) plane that 
has more than one count at the exact same real-valued po- 
sition (which would be the case if a pixel contained more 



Chandra Interactive Analysis and Observations; developed and 
maintained bv the Chandr a X-ray Center, and available at 



references therein) implemented in the Chandra x-ray data http://asc.harvard.edu/ciao 
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than one count) becomes disproportionately significant, and 
distorts the value of Dks- 

Therefore, to apply a 2D K-S test to 3D images, it is first 
necessary to "un-bin" the images by spreading the counts in 
each pixel over the area of the pixel. Since no information 
is preserved about the exact impact locations of the pho- 
tons on the detector, we introduce a random shift (between 
—0.5 and +0.5 pix) in the integer- valued coordinates of each 
count. Every count in the image is thus assigned a unique po- 
sition (within the precision limits of the computer) , and the 
integer-valued 3D coordinates are converted to real-valued 
2D coordinates. Since this routine distributes the data on a 
sub-pixel scale, we refer to it as "subpixelization." 

Unfortunately, subpixelization incurs an undesirable ef- 
fect, due to the randomness with which the counts are moved 
around within the pixels. In particular, two subpixelized ver- 
sions of the same image are never the same. Therefore, K-S 
testing of different subpixelizations of the same two images 
will produce different results for Zn every time. Our expe- 
rience is that Zn varies by about 2 to 4 per cent between 
runs; the corresponding variations in the significance level 
SL may be as high as ±8 per cent for values of Zn near 
SL = 50 per cent (Table To obtain a mean value for the 
K-S probability P(> Zn) (equal to 1 — SL) with which two 
images represent the same parent distribution, it is neces- 
sary to run the 2D K-S test multiple times. The error in the 
K-S probability can be estimated from the variations in Zn 
among the different runs. 

The above method may appear contrived and unneces- 
sary, when instead of applying a 2D K-S test, by following 
the generalisation in FF, a 3D K-S test can be implemented. 
A 3D test does not have the undesirable uncertainties associ- 
ated with subpixelization, and gives the exact K-S probabil- 
ity that two images represent the same parent distribution. 
FF report that the 3D test exhibits greater power in reject- 
ing wrong hypotheses than a three-fold 2D test along each of 
the {x,y), {y,z) and {z,x) planes. However, our experiments 
with the tests on simulated ROSAT images show that the 
2D K-S test applied to subpixelized images is more powerful 
than its 3D counterpart (with accordingly generated look-up 
tables) applied to the original (not subpixelized) images. K- 
S tests performed on simulations of an 880-count five-source 
distribution with centre-to-centre distances of 1.1-2.1 times 
the PSF size (5 arcsec for ROSAT/RRl) show that the 2D 
test on subpixelized images can distinguish individual source 
position shifts as small as 2-3 arcsec (depending on the rel- 
ative source locations and intensity) at the > 97 per cent 
significance level. The 3D test finds such shifts insignificant: 
it gives a probability of > 70 per cent that the simulations 
represent the same parent distribution. 

The higher power of the subpixelized 2D test is ex- 
plained by the manner in which pixels are weighted. Subpix- 
elization ensures that all photons are assigned equal weights: 
a desirable effect, since pixels are weighted proportionally to 
their intensity. The 3D K-S test, on the other hand, assigns 
equal weights to all pixels, so single-count pixels (often from 
background) are as significant as pixels with multiple counts 
(denoting sources). Hence, the 2D test on subpixelized im- 
ages is more sensitive to variations in the source distribution 
than its 3D counterpart. 

We also note that subpixelized images represent more 
realistically the incoming photon distribution. Repeated 2D 



K-S tests on independent subpixelizations produce an esti- 
mate of the random error in the K-S probability induced by 
photon binning. 



3.1.2 One-sample vs. two-sample two-dimensional K-S 
test 

When comparing an image to a proposed model distribu- 
tion of sources (using an analytical or a fitted PSF), perfect 
information is available about the shape of the model dis- 
tribution. It is therefore appropriate to use a one-sample 
K-S test to compare the image to the model distribution. 
However, the analytical form of a fitted PSF is often very 
complex. The added complication of having several nearby 
sources with overlapping PSFs (as in a crowded field) makes 
it computationally very tedious to calculate the fraction of 
the analytic model distribution in the quadrants around ev- 
ery count i n th e image (needed to compute the Dks statis- 
tic; Section 



2.1) 



To avoid lengthy 2D surface integrals, the two-sample 
K-S test can be used instead, to compare the image to a 
simulation generated from the proposed model. The ran- 
dom deviations in the representation of the model can be 
decreased if a "bright" (high number of counts n) simula- 
tion is used, that follows the model closely. The fractional 
deviation of the simulated vs. expected counts per pixel will 
decrease as l/-yncts for large nets, where ricts is the number 
of counts per pixel in the model image. The running time 
of the two-sample 2D K-S test is however proportional to 
the square of the sum ni -I- 712 of counts in the samples com- 
pared, and it is impractical to use the two-sample K-S test 
on simulations containing high number of counts. 

To take advantage of the higher power of the one- 
sample test (given perfect information about the shape of 
the model distribution), and to avoid long running time, we 
compare the image (containing n counts) to a bright sim- 
ulation (njnod = k X n counts; k ~ 50) of the proposed 
model, using the one-sample test. In essence, we use Monte 
Carlo integration for the model, the assumption being that a 
bright simulation can be made to represent the model with 
sufficient accuracy, and simple summing of the counts in 
the quadrants can be substituted for analytic integration of 
the PSF. We therefore do not expect the properties of this 
pseudo one-sample K-S test to be significantly different from 
those of the FF one-sample K-S test. In particular, we as- 
sume that the new test is still distribution-free for a fixed 
correlation coefficient CC of the bright simulation (hereafter 
referred to as the "model"), and that its properties vary 
slowly with CC (as observed in FF for their 2D K-S test). 
Clearly, our test will converge to the FF one-sample test for 
k ^ 00, but will run faster than the latter for moderate- 
sized k, since lengthy analytical integrations are substituted 
with Monte Carlo integration. The pseudo one-sample test 
(running time proportional to n{n -\- nmod) = (A; -I- l)n^) is 
also fe -I- 1 times faster than the two-sample test (running 
time proportional to (n -I- rimod)^ — {k -\- l)^n^). Thus, for 
moderate-sized k, the pseudo test approximates the power 
of the one-sample K-S test, and is faster than both the one- 
sample and the two-sample K-S tests for a general model 
distribution. 

Table ^ presents a comparison of the performance of 
the two versions of the test, using the same five-source 
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Table 2. Comparison of the performance of the pseudo one-sample 2D K-S test and the 



FF two-sample K-S test as a function of fc = ^ 


for a model with CC = 0.10. 






pseudo one-sample test 






two-sample test 




k 


ni 




P(> Zn) 






P(> z^Y 


1 


880 


1.78 ±0.07 




440 


1.30 ±0.04 


34 ± 5% 


2 


880 


1.72 ±0.05 




587 


1.44 ± 0.03 


19 ± 3% 


5 


880 


1.56 ±0.04 




733 


1.48 ± 0.04 


16 ± 4% 


10 


880 


1.88 ±0.05 




800 


1.79 ±0.03 


3±1% 


20 


880 


1.59 ±0.06 




838 


1.64 ±0.03 


8 ±2% 


50 


880 


1.62 ±0.04 


<2%'* 


863 


1.65 ±0.05 


8 ±3% 


100 


880 


1.62 ±0.06 




871 


1.68"= 


~7% 


200 


880 


1.66 ±0.06 




876 


1.74= 


~ 6% 



° Errors determined from K-S comparisons of 10 different subpixelizations of the same 
simulations. Systematic uncertainties associated with representing a continuous model 
by a discrete distribution are not included. 

^ For the two-sample test n = "^"^ = — — = -r^-r'n.l- 

" Obtained from Table |l] for CC = 0.10 0.0 of the two samples. 

Interpolated from lines 1 (fe » 50) and 2 (fc = 25) of Table | (CC = 0.10 Ri 0.12). 

Only one comparison was performed due to the longer running time. 



Table 3. Critical values Z„_sL for simulated source distributions. Model size rtmod = 45000 ~ 
50 X 880, i.e. k » 50. 







5L(%)t 


30 


40 


50 


60 


70 


80 


90 


95 


99 


CC 


n* 


# of simul. 




















0.12 


880 


10000 


0.98 


1.03 


1.08 


1.13 


1.20 


1.27 


1.38 


1.47 


1.67 


0.12 


1800 


10000 


1.03 


1.08 


1.13 


1.19 


1.24 


1.32 


1.43 


1.54 


1.74 


0.23 


820 


10000 


0.95 


1.00 


1.05 


1.09 


1.15 


1.22 


1.33 


1.43 


1.61 


0.23 


1800 


10000 


1.02 


1.07 


1.11 


1.17 


1.23 


1.30 


1.42 


1.52 


1.71 



t Significance level {SL = 1 - P(> Z„)). 
t Size of sample. 



setup as in Section with one source shifted by only 

0.4 FWHMs = 4 pix between the model and the simula- 
tions. As expected, the power of the two-sample test in- 
creases with increasing fc, but for fc > 20 the results for 
Zn,2a do not chauge significantly. The values of Z„,is from 
the pseudo one-sample test exhibit greater random varia- 
tions for small fc than those of Zn,2s- These are due to the 
fact that the small-number statistical uncertainties in the 
two samples are not averaged out in the one-sample test, 
whereas in the two-sample test they are. For fc ^ 20 how- 
ever, the values for Z„,is become self-consistent (as well as 
consistent with those from the two-sample test), and there 
is little power to be gained from increasing fc. 

Despite the similar behaviour of the FF K-S test and 
our pseudo one-sample test, the reference tables for the for- 
mer (Tables A1-A5 in FF; Table ^ cannot be used for the 
latter, since the Z„ distributions are different. We there- 
fore ran Monte Carlo simulations to create a separate look- 
up table for the new test. Table ^ is based on comparing 
simulated source distributions with ni ~ 850 counts to 
fc « 50 times brighter simulations ("models," n2 ~ 45000) 
of the same source distributions. Sources were simulated on 
a 128 X 128 pix field using the analytic 5-arcsec (10-pix) HRI 
PSF. We also ran Monte Carlo simulations with higher ni 
(ni ~ 1800) against « 45000-count models for the pur- 
pose of interpolation. Although in this case the value of fc 
is lower (fc = 25), we choose to set the value of n2 = 45000 
as the standard, as it is a measure of the precision of the 
Monte Carlo integration. High values of the sample corre- 



lation coefficient CC were not pursued, since astronomical 
images rarely have highly correlated photon distributions 
(esp. in crowded fields, but also given random background). 
Thus constructed. Table ^ should be applicable to most low 
background, low S/N imaging cases (esp. x-ray). 

Using Table ^ for the values of Zn^sL, the results for the 
K-S probability (P > Zn) from the pseudo one-sample K-S 
test will be consistent with those from the FF K-S test using 
the FF tables. Hereafter, unless explicitly stated otherwise, 
we shall refer to our pseudo one-sample 2D K-S test as "the 
2D K-S test," or simply as "the K-S test." 



3.2 Application of the two-dimensional K-S test 
to simulated images 

After establishing that the 2D K-S test can be successfully 
applied to images, it is important to determine the sensi- 
tivity of the test with respect to deviations from the pro- 
posed model. Here we test the responsiveness of the 2D K-S 
test with respect to changes in the parameters of individual 
sources in the model distribution. We look for the minimum 
deviation in a single parameter (position or intensity of a 
source) that enables the test to tell the distributions apart 
at the > 99 per cent significance level. As a trial source 
distribution we use the one simulated in Figure whose 
parameters are listed in Table ^ (sources X1-X5), and com- 
pare it to models that have one parameter changed. 

It is worth noting here, that given the significant over- 
lap of PSFs in a crowded field such as the one in Figure Sa, 
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Table 4. Parameters of the source distribution best describing 
the x-ray emission in the central 128 X 128 pix section of the 
NGC 6397 ROSAT image. 



Table 5. Number of sources in the best-fit model vs. K-S prob- 
ability. 



Source 

XI 
X2 
X3 
X4 
X5 
X6 
bkg 



Pixel 



coordinates 



(x - 4012) (y - 4027) 



71 

56 
59 
73 
76 
94 



55 
56 
37 
39 
88 
77 



Counts 

78 
176 
161 
88 
197 
16 

0.015 cts/pix 



Counts 
(x50) 
3900 
8800 
8050 
4400 
9850 
800 
0.75 cts/pix 



1 (a) : 

xa XI 

X3 X4 


-(b) . : 

• ; 


1 1' 1 1 1 1 1 1 1 1 1 1 

- (c) ^ : 


1 ' 1 1 1 1 1 ; 1 > ' 1 

- (d) ^ : 



50 

X - 4012 



50 100 
X - 4012 



Figure 2. Scenarios for determining the sensitivity of the pseudo 
one-sample K-S test with respect to changes in the parent dis- 
tribution. The X- and jz-axes are HRI pixels. A 10-pixel FWHM 
PSF has been used, (a) A subpixelized 880-count 5-source simu- 
lation. The source parameters are listed in Table ^ (X1-X5). The 
simulations in (b)-(d) are 50 times brighter with a single source 
modified in each: (b) source X2 is moved 2 pix up and 2 pix to the 
left; (c) source X2 is 1.45 times brighter (11600 cts); (d) source 
X2 is at 60 per cent intensity (4800 cts). 



the sensitivity of the test with respect to changes in the 
distribution need not be isotropic. It depends on the rela- 
tive positioning and brightness of the sources. We therefore 
test several possible scenarios of such changes. The scenarios 
shown in Figure |^ and described in its caption correspond 
to the limiting cases, in which the K-S test can distinguish 
the distributions at the > 99 per cent level. 

After performing K-S tests between the simulation in 
Figure ^a and the models in Figures ^-d, we establish that 
the positional sensitivity of the K-S test in crowded fields 
depends on the relative source intensity, and on the direc- 
tion in which a source is allowed to move with respect to the 
crowded region. The sensitivity to moving a source is higher 
for brighter sources (3 to 4 pix for source X2) than for the 
fainter ones (4 to 9 pix for source XI), and is generally (al- 
though not conclusively) poorer when the source is moved 
toward the region of crowding (9 pix for source XI, 3 pix 
for X2) as opposed to when it is moved away (4 pix for XI 
and X2). The brightness sensitivity of the K-S test is also 



N 
3 
4 
5 
6 



CC of model 
0.21 
0.22 
0.19 
0.21 



1.88 ± 0.04 
1.48 ± 0.05 
1.23 ± 0.05 
1.19 ± 0.04 



K-S prob P(> Z„)" 
< 1% 
4±1% 
22 ± 6% 
26 ± 6% 



p/ b 

1.5% 
43% 
90%^= 
(100%) 



" Forn = 980 counts in the NGC 6397 image. 
^ Probability that a best-fit model with TV sources and the 
ROSAT image represent the same parent distribution; 1 — 
is the significance of a ddi ng an (A' -|- l)-st source. P'j^ is calcu- 
lated as Ppf 5 (Section 7.2). 

Estimated as Ps^e, i.e. the probability that a 5-source best-fit 
distribution can represent a 6-source one; 1 — P5.6 = 10 per cent 
is the significance of adding a sixth source. 

dependent on the relative source intensity in a crowded field: 
brighter sources can vary by a smaller percentage (~ 45 per 
cent for X2) than fainter sources (~ 70 per cent for XI). 
For a source that is sufficiently far away (~ 3 FWHM) 
from the crowded region (source X5) the sensitivity of the 
K-S test is greater and nearly position- independent. For 
source X5 we set the limits at a 2 pix shift in any direc- 
tion (P(> Zn) ~ 1 per cent) or a 30 per cent change in 
intensity (P(> Zn) < 1 per cent). 

The above-determined sensitivity limits are based on 
varying only one parameter (position or intensity) of a sin- 
gle source, while keeping all other parameters fixed. This is 
the approach used to determine the 95 per cent (~ 2a) con- 
fidence limits (Table ^) on the positions and intensities of 
the detected sources in our NGC 6397 image (Section H). 



4 THE 2D K-S TEST IN PARAMETER POINT 
ESTIMATION FOR SOURCE DETECTION 

4.1 Algorithm 

An iterative source-fitting algorithm was devised that aims 
to minimise the Zn statistic, thus maximising the proba- 
bility that an image and a simulation represent the same 
parent distribution of sources. The final simulation that re- 
sults from this algorithm will contain the best estimate for 
the number, positions and intensities of the sources in the 
image, subject to limitations arising from the sensitivity of 
the test. The iterative procedure steps through the following 
algorithm: 

(i) An initial guess of the source distribution is made. 
This can be a source at the location of the brightest pixel 
(thus starting with a one-source configuration) or a guess 
with A'initiai > 1 number of sources. Both initial guesses will 
produce the same results for a distribution with A^flnai 5S 
initial sourccs. The PSF is fit to a single unresolved source 
in an uncrowded part of the image so that aspect or other 
image systematics are included. 

(ii ) A bright simulation (a.k.a. a "model"; see Sec- 
tion 3.1.2) based on the current guess for the source dis- 
tribution is created (using the fitted PSF) and normalised 
to the image intensity. The normalised model and the image 
are smoothed with a Gaussian function to roughly match 
the ROSAT/HRl resolution. The residual between the two 
is then formed. 
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Figure 3. Values of Z„ vs. the maximum deviation from zero, 
f max, in the smoothed residual. The error-bars in Z„ represent 
the standard deviation of Z„ due to 10 independent subpixeliza- 
tions of the compared simulations. The statistics Zn and Dmax 
are highly correlated (correlation coefficient 0.93), and we there- 
fore use a minimum in i)max as an indication of being near a 
minimum in Zn (and hence, for constant CC , near a maximum 
in P(> Z„)). 



(iii) The next guess for the source distribution is obtained 
by moving the source positions against the steepest gradient 
in the residual, and by adjusting the intensities, so as to 
decrease the maximum deviation from zero (Dmax) in the 
smoothed residual. We indeed observe that Dmax is strongly 
correlated to the value of Zn (correlation is 0.93; Figure 3f), 
and hence (for constant CC) to the K-S probability P{> 

Zn). 

(iv) Repeat steps (ii) and (iii) until Dmax is minimised. 
Compare the image to the model with the 2D K-S test and, 
if necessary, further minimise Zn by applying small changes 
(e.g. single-pixel shifts) to the model (since the nature of the 
relation between Dmax and Z„ is not established rigorously) . 
The final simulation will contain the best guess of the po- 
sitions and intensities of the assumed sources. The image is 
compared to the model using the 2D K-S test. 

(v) Steps (i) through (iv) are run for a fixed number of 
sources (guessed in step (i)). If the final K-S probability of 
similarity is not satisfactory (e.g. not > 5 per cent), a new 
source is added at the location of the highest residual, and 
the algorithm is repeated from step (ii). 

(vi) If the addition of the last source did not incur a de- 
crease in the Zn statistic larger than its uncertainty (±2 per 
cent to ±4 per cent), the last added source is considered 
marginal, and the previous best guess for the number, posi- 
tions and intensities is taken as the final one. 



Figure 4. A comparison between a Zt&, a^nd a Zt^„ curve ob 



As noted in step (i) an initial guess with TVinitiai > 1 
number of sources can also be fed into the algorithm. Such 
a guess can be made either from visual inspection of the im- 
age, or after applying a deconvolution algorithm. We found 
that L-R deconvolution (see Section ^ below) gives good ini- 
tial estimates of the positions of the individual sources. How- 
ever, since deconvolution can introduce spurious sources, the 
initial guess of the number of sources should be conservative. 



880 

tained using the 5-arcsec predicted HRI PSF. The ^ggQ curve is 
obtained from K-S tests between a 5-source model and 10000 real- 
isations of a 4-source simulation that best represents the 5-source 
model. The overlap area P4 5 determines the false detection prob- 
ability of the 5th source. Here P4,5 = 0.23. 



4.2 Performance 

The above procedure has not been automated, and therefore 
(due to subjectivity in "guessing" a simulated source distri- 
bution after having created it) we have not performed tests 
to explicitly determine its efficiency in detecting sources. We 
quote the ability of the 2D K-S test to detect small changes 
in the positions (within ~ 0.2 FW HM) and/or intensities 
of individual sources (Section 3.2) as an indication of the 
power of the iterative algorithm. Nevertheless, we have de- 
vised a method to test the confidence with which a certain 
number of sources can be claimed in a given photon dis- 
tribution. The method takes our best guess for the source 
distribution in the image with a given number of sources 
(e.g. A'^), and compares a model of it to a faint simulation of 
our best guess with one source fewer (TV — 1). In this way we 
can test in what fraction Pn-i,n of the cases our proposed 
model (with N sources) can describe a source distribution 
with N — 1 sources. In other words, we test for the signifi- 
cance (1 — Pjv_i^jv) of the addition the N-ih source; Pn-i,n 
is thus its false-detection probability. If this comparison is 
performed many of times (of the order of the number of 
Monte Carlo simulations done for each row in Table a 
^N-i,N g^jj-yg fQj. ^jjg guesses is recovered. The latter 
can be then compared to a Zn'^ curve, obtained in a simi- 
lar fashion comparing A-source simulations to an A-source 
model. 

Example Zn'^ and Zn~^'^ curves for N — 5 and n — 
880 are shown in Figure ^ The overlap of the two curves 
gives the fraction of Monte Carlo simulations, in which a 
best-fit distribution with A — 1 sources produces an image 
that has the same K-S similarity to the model as that of a 
best-fit distribution with A sources. The ratio of the overlap 
area to the area of either of the Zn curves (assuming they 
are both normalised to the same area) is the desired fraction 

f JV-1,JV. 

To investigate the dependence of the overlap area 
Pn~i,n on the width of the PSF, we ran K-S tests simula- 
tions built with an 8-arcsec Gaussian PSF. The result is that 



for a wider PSF the Zf! 



curve is narrower and is shifted 
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toward smaller Z„ (the model and the simulations look more 
alike). The Z^''^ curve however is not affected. The overall 
effect is that the false-detection probability Pn-i,n of the 
A''-th source increases (from 0.23 to 0.43 for N = 5; Table ||). 

The above technique can be generalised to produce 
fjv-i,jv for arbitrary integers i < N and j. An application 
of Pjv-i,]v is discussed in Section 7.2 



APPLICATION TO A DEEP ROSAT IMAGE 
OF NGC 6397 



The iterative source-modelling procedure (Section [4.l[ ) was 
applied to our 75 ksec ROSAT/URl exposure (March 1995) 
of the core region of the post core-collapse globular cluster 
NGC 6397 (Figure 3). Standard asp ect correction routines 



(Harris et al. 1998a, Harris 1999a b[) were applied to the 



image to improve the S/N ratio. After the aspect corrections 
the PSF improved from 10.3 arcsec x 8.3 arcsec to 8.3 arcsec 
X 7.9 arcsec, as measured from the shape of a background 
point-source quasi -stellar object ( QSO) at 3.7 arcmin off- 
axis (source "D" in Cool et al. 1993). The obtained size of the 
PSF was still much worse than the predicted 5 arcsec. This 
effect is not due to the known deterioration of the ROSAT 
PSF with increasing off-axis angle beyond ~ 5 arcmin, since 
the QSO is only 3.7 arcmin from the centre of the field. 
Residual (unknown) aspect errors are present, and in the 
analvsi s below we use a fitted PSF instead of the nominal 



We analyse the central 128 pixel (64-arcsec) square re- 
gion of the image, containing 980 counts. Model simulations 
were created using an analytical PSF fit to the QSO with 
the iraf/daophot routines PSF and ADDSTAR. The PSF 
was comprised of a FWHM ~ 8-arcsec Gaussian core and 
Lorentzian wings, where the core and the wings could be 
tilted along different directions in the image. In determining 
the false-detection probabilities (from the overlap of the Zn 
curves) however, for faster iteration we used a symmetri- 
cal 8-arcsec Gaussian PSF, noting that the Zn distributions 
based on the fitted PSF and on the Gaussian PSF are ex- 
pected to be indistinguishable, since the test is distribution- 
free. 

Table ^ presents results for the K-S statistics of the 
best-fit models for a given number of sources. The error in 
the val ues of Zn and P(> Zn) is the one-sigma uncertainty 



4, 5 and 6 sources are shown in Figure |^. The number of 
counts per source for the five- and six-source cases are the 
same as listed in Table ^ 

Source X6 is sufficiently faint and detached from the 
central group (X1-X4), that its addition did not necessi- 
tate any changes in the prior five-source configuration. It is 
11 arcsec away from the closest source (XI), and 4.5 times 
dimmer than the faintest one (also XI; Table |). The false 
detection probability (determined as Ps.e) for source X6 is 
90 per cent. 

Source Xi,4 in the four-source best-fit solution falls ap- 
proximately in the middle between sources XI and X4 of 
the five-source solution, which is consistent with them hav- 
ing conxparable intensities and being fainter than X2 and X3 
(Table g). 

The derived optimal number of five sources in our 
75 ksec im age of NGC 639 7 is consistent with the earlier 
suggestion (Cool et al. 1993) that at least four x-ray sources 
(XI, X2, X3 and X5) are present in an 18 ksec exposure of 
the same region (found by visual inspection of the peaks in 
the image, and confirmed by a one-dimensional azimuthal 
K-S test around the source centroids). The locations of the 
detected sources are also consistent (up to a ~ 2-arcsec 
systematic offset) with the positions of known optical cata- 
clysmic variables (CVs): the optimum five-source positions 
(solid triangles. Figure ^ are systematically displaced by 
~ 2 arcsec to the lower right of CV1-CV5 (open t riangles, 
as measured from Ha, TZ, and UBVI HST data; Cool et 



al. 1995, 1998). Such an offset is well within the expected 



variation (< 5 arcsec) of the absolute pointing of ROSAT. 

Source X5 does not have a known optical counter- 
part, and sources CVl and CV4 are unresolved in the 
ROSAT x-ray image due to their proximity (~ 2.5 arcsec 
~ 0.3 FWHM), so that the emission of XI is probably due 
t o both of them . The object CV5 was previously known 
(| Cool et al. 1998| ) as a near t/V-excess star, but proposed as 



a CV-candidate (Crindlay 1998; Metchev 1999) only after 



re-evaluating its probable association with the x-ray emis- 
sion from X4. Its subsequent confirmation as a CV-candidate 
( Hg-emission objec t) in our foUowup deep HST imaging 
(Taylor et al. 2001) is indicative of the reliability of our 



iteration scheme to detect faint sources in crowded fields. 

Similarly, source X6 (if real) may be assoc iated with 
the Chandra source U43 (Grindlay et al. 2001c): a proba- 



ble faint BY Draconis binary, identified as PC-4 in Taylor 



due to f ubpixelization, as determined from K-S comparisons et al. (2001). Such objects are expected in globulars as the 



between the model and 10 independent subpixelizations of 
the same image. It is an estimate of the error in the mean of 
the Zn distribution of comparisons between the image and 
the A'^-source model. 

Following the logic of step 6. in the K-S probability 
maximisation algorithm (Section |4.l[ ) , we conclude that four 
sources are insufficient to represent the image conclusivelv. 



binary companions of stars on to which they have trans- 
ferred their envelopes , and indeed the reported velocity in 
Edmonds et al. (1999) suggests a massive but dark binary 



companion. The latter is most likely a neutron star (NS), 
since He-WD/NS systems are expected in millisecond pul- 
sars (MSPs), and MSPs are detected as faint x-r ay sources 
with luminosit ies comparable to that of X6 (cf. Becker & 



since tl e addition of a fifth source decreases the Zn statis- 



tic by more (17 per cent) than the 3.5 per cent error in 
the Zn statistic for four sources. However, the addition of 
a sixth source is not justified, since the decrease (0.04) in 
Zn is smaller than the error (0.05). We therefore claim that 
five sources are sufficient, and that at least four sources are 
necessary (at the 1 — Pa, 5 > 99 per cent level) to account for 
the observed photon distribution in the ROSAT ima.ge. The 
source centroids for the optimal source configurations with 



Triimper 1999| ). 

Table |6| lists results for the detected sources. A mean 
bore-sight offset of — 3.5±1.0 pix in x, and -|-1.7±1.0 pix in y 
(-1.8 arcsec and +0.85 arcsec, respectively) has been applied 
to the x-ray positions to align them with the suggested opti- 
cal counterparts. In calculating the offset, we discard source 
XI (since the x-ray emission in its vicinity comes most likely 
from both CVl and CV4), and weigh the measured shifts in- 
versely to the square of the positional uncertainties of the 
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Figure 5. A 75 ksec ROSAT/^iRl exposure of the central re- 
gion of NGC 6397, with the detected sources marked. The x-ray 
image is smoothed with a 2-d a = 2 arcsec Gaussian. The clus- 



ter cent -e is at (etc, 5c) = (17:40:41.3, -53:40:25), (Djorgovski & 



Meylan 1993). The conversion from pixel to celestial coordinates 



is accurate to within 1 arcsec. 

sources along x and y. The 95 per cent confidence radius of 
the source coordinates corresponds to the maximum shift of 
a single source from its listed position (keeping all intensi- 
ties and other source positions constant) that maintains the 
K-S probability above 5 per cent. 



6 COMPARISON TO OTHER SOURCE 
DETECTION METHODS 

To get an idea of the superior performance of our source- 
modelling scheme we compared it to established source- 
detection algorithms, such as the classical "sliding-cell" de- 
tect, the wavelet detect, the iraf/daophot PSF-fitting task 
ALLSTAR, and ML analysis. We also used deconvolution rou- 
tines on the image to determine possible source locations. 
Below we discuss briefiy each of these alternatives. 

Sliding-cell Detect: The sliding-cell detect algorithm is 
based on S/N calculation and was not expected to perform 
well in a crowded low-5'/A*' field. Indeed, the two versions of 
this algorithm in the iraf/pros package (tasks imdetect 
and LDETECT in the xspatial package) fail to produce the 
expected number of x-ray sources in the cluster, imdetect 
uses a constant average background for the entire image, and 
a variable detect cell size (squares with sides from 4 arcsec 
to 24 arcsec) to search for sources. The larger detect cells 
fail to find more than three sources in the central region of 
NGC 6397, whereas the 4-arcsec detect cell size is too small 
for use with our PSF (FWHM ~ 8 arcsec), and produces 
an unjustified high number of individual sources. The lo- 
cal det ect algorithm (task ldetect) calculates SIN around 



Figure 6. Best-fit source positions for 4, 5 and 6 sources. Open 
triangles mark the positions of known H g-emission object s (prob 



able CVs. numbered with their ID from Zoo\ et al. 1995: Urind 



lay 199£ ) - candidate c ounterparts for s ources X1-X4. Open 



squares mark (71/-excess (Z^ool et al. 1995; Edmonds et al. 1999) 
and/or f aint (L^, < 10'^^ erg s~^) Chandra sources (IDs from 
Fig. 1 of [Srindlay et al. 2001a ). The consistent shift between the 
HST/ Chandra CV and the RO S AT :k-i ay positions is indicative 
of a bore-sight offset (~ 2 arcsec) of ROSAT for this observation. 



1.5 arcsec and 2.5 arcsec from the source) as an estimate of 
the noise. As a result it does not handle crowded fields ad- 
equately, and cannot distinguish blended sources. Even the 
two smallest detection cells (6 arcsec x 6 arcsec and 9 arcsec 
X 9 arcsec) do not find more than 3 sources in the image in 
Figure ^ 

Wavelet Detect: This algorithm based on the wavelet 
transform has only recently b e en applied to imagin g as- 



tronomy (Freeman et al. 2001 



Damiani et al. 1997, and 



references therein), and has been demonstrated to outper- 
form other source detection algorithms in \ow-S/N fields. 
We used an implementation of the wavelet detect based on 
the Marr wavelet, or the "Mexican Hat" function, coded in 
the wavdetect task in the Chandra detect 1.0 Package. 
The algorithm is most sensitive to structures of size approx- 
imately equal to the width of the Mexican Hat function. 
Running wavdetect on our NGC 6397 image (Figure ||) 
with transforms of width 8 arcsec produced only the same 
three sources already found by the sliding-cell algorithms. 
This was not unexpected, since in simulated images for the 
Chandra High Resolution Camera (FWHM — 0.5 arcsec), 
wavdetect is unable to discern point sources less than 2 
FWHM apart.| 

Image Deconvolution: There exist a number of widely 
used image deconvolution algorithms that are applicable to 



Chandra DETECT User's Guide: TTRT.: |^ttp://hea 



each pi: eel, using the local background (in a region between 



www.harvard.edu/asclocal/user / swdocs / detect /html / 
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Table 6. X-ray sources detected in NGC 6397. 



Source 


R.A. 


DEC. 


95% conf. 


Count rate 


(erg/s)« 


Optical 


Chandra offset** 


(arcsec 




(J2000) 


(J2000) 


radius 


(ksec^^ ) 


kT = 10 keV 


counterpart^ 


ajs — ac - 


XI 


17:40:41.46 


-53:40:18.4 


2.5" 


1.2 ± 0.5 


3.3 X 10^1 


CVl, CV4 


-0.2 


-0.6 


X2 


17:40:42.48 


-53:40:18.0 


1.5" 


2.7 ±0.6 


7.4 X lO^i 


CV3 


0.7 


1.2 


X3 


17:40:42.32 


-53:40:27.4 


1.5" X 2.5" 


2.5 ± 0.6 


6.9 X lO^i 


U18, CV2 


0.6 


-1.8 


X4 


17:40:41.53 


-53:40:26.4 


2.0" 


1.3 ± 0.5 


3.6 X lO^i 


CV5 


0.7 


-0.1 


X5 


17:40:41.36 


-53:40:02.0 


1.2" 


3.0 ±0.6 


8.2 X 10^1 


U24 


0.0 


0.0 


X6 


17:40:40.4 


-53:40:18 


oo'' 


0.4 


< 1.1 X 10^1 


U43 


0.8 


2.8 


" Unabsorbed luminosities listed in 


the 0.5-2.5 keV band, for a 


cluster distance of 2.2 koc. and column density of 1.0 



lO^Lcm ^. For best fit {Chandra) column densities and bremsstrahlung spectra for individual sources, see Grindlay et al 



(2(|l01c|). 

^ The 2(7 confidence radius of the position of X6 is infinite, because the K-S probability that the model and the image represent 
the same parent distribu tion is always above 5 per cent, regardless of the source location. 



U18. U2 4 and U43 are [Grindlay et al. (2001^ ) Chand ra IDs. U18 also ide ntified as either a BY Dra or MSP by [Grindlav et 
al. | (2001e ), and U43 identified as a BY Dra binary by Taylor et al. (2001). CV2 first identified as Ho object by |Cool et al, 

(i<g^ 

Offset between the given positions (subscript R) and the ones listed in Grindlay et al. (200l£ , subscript C). A boresight 
offset {Aaji_c = 4.9 pix and ASji_c = 4.6 pix) has been applied to match the positions of our best-constrained source 
(X5) and its Chandra counterpart (U24). When the emission from two Chandra/HST sources corresponds to a single ROSAT 
source, the latter has been associated with the mean position of the Chandra/HST sources. 



moderately crowded fields. After comparing results from the 
IRAF implementations of the Maximum Entropy Method, 
the L-R algorithm (both applicable primarily to optical im- 
ages), and from CLEAN (used mostly in radio imaging), we 
found that L-R deconvolution (Lucy 1974; Richardson 1972) 



most reliably discerns the five-source distribution found by 
our iterative source-modelling scheme (Section ^). The po- 
sitions of the peaks in the deconvolved image are in excel- 
lent agreement (to within ±1 pix = ±0.5 arcsec) with the 
K-S best-fit source positions, which exemplifies the useful- 
ness of L-R deconvolution in analysing crowded fields. Un- 
fortunately, the L-R method does not provide a measure 
of the goodness of fit of these positions and of the signifi- 
cance of the peaks in the reconstructed image. These need 
to be determined separately with a multi-source fitting rou- 
tine (since the field is crowded), such as daophot/allstar, 
or the current (2D K-S) iterative method. Furthermore, the 
obtained intensities of the deconvolved sources are in much 
poorer agreement with the ones from the 2D K-S best-fit 
model. Nevertheless, L-R deconvolution does give an indica- 
tion for the existence of more than 3 sources (which could 
not be determined with the source-searching methods) . The 
L-R method thus provides a very good initial guess for the 
source configuration, which can be input to iterative source- 
modelling algorithms. 

daophot/allstar: The daophot package is designed 
for the analysis of crowded optical images, and as such 
it assumes that the images are in the Gaussian statistics 
(high number of counts per pixel) regime. Thus, strictly 
speaking, the package is inapplicable to data governed by 
Poisson statistics, such as most x-ray images (including our 
NGC 6397 image, containing ^ 3 counts per pixel), because 
it severely underestimates random errors. However, until re- 
cently DAOPHOT was the only widely available softwar e for 
reduction of crowded fields, and it has been suggested (dool 
et al. 1! 93) that it can be useful for analysing crowded x-ray 
fields. 

Our experience with ALLSTAR is that it is heavily depen- 
dent on several loosely defined parameters which, in regimes 
of severe source confusion and low signal-to-noise as in our 



NGC 6397 image (Figure g), critically determine the per- 
formance of the task. We found that different combinations 
of the values of the parameters and of the initial guess for 
the source distribution produced different final results, in 
which the number of detected sources in the NGC 6397 im- 
age varied from 2 to 5. By judiciously adjusting its parame- 
ters, ALLSTAR can be made to detect 5 sources, however that 
combination is not favoured statistically over other combina- 
tions with fewer sources. In the case when ALLSTAR detects 5 
sources, the obtained positions and intensities are such that 
the K-S probability of similarity with the ROSAT image is 
< 1 per cent {Z„ = 2.2). 

Maximum Likelihood: Given our method of optimisation 
- minimising the maximum residual Dmax (albeit we then 
further minimise the K-S statistic Zn) - ML analysis would 
be expect ed to produce a similar fit. This is indeed the ap- 
proach of Verbunt & Johnston (200C ) in analysing the same 
RO SAT field. The rcsuhs for t he 5 detected sources (Model I 
in (Verbunt & Johnston 2000); X1-X5 in this paper) agree 
well; in addition, our analysis suggests the possible presence 
of the faint source X6. We choose to employ a 2D K-S test to 
assess the goodness of fit instead, banking on its sensitivity 
to diffuse distributions. As pointed out by the referee, it is a 
good test for the location of smeared objects, but it is rather 
insensitive to their width. Via K-S, a source may be deduced 
to be unresolved, despite having broader profile, which can 
frequently be the case in Poisson noise limited images. 

We have thus demonstrated that under conditions of 
severe source confusion and low S/N, our source-detection 
method based on a 2D K-S test works better than other 
available techniques. We attribute its performance to the 
fact that our approach uses the actual PSF in searching for 
sources (sliding-cell and wavelet detect algorithms do not), 
that no information is lost to binning (as in the Pearson 
test, used in ALLSTAr), and that it is more sensitive to broad 
emission than other tests (e.g., ML). 

We have not made a compariso n of our method agains t 
the Pixon deconvolution method (Pina & Puetter 1993). 
Our method was originally intended to enhance sensitivity 
for crowded point source detection; the Pixon method also 
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shows good results for the detection of low surface brightness 
features. 



7 DISCUSSION 

7.1 Applicability to distributions with unknown 
parameters 

Rigorously, the presented look-up Table ^ (generated by 
comparing simulations of models with a priori known pa- 
rameters) is not applicable when comparing an image of an 
unknown source distribut ion to a simulation with known pa- 
rameters. 



Lilliefors (1967) investigates this situation for the 



case of the ID K-S test and sampling from a distribution 
with unknown mean and variance ("the Lilliefors test for 
Normality"). He finds that the standard ID K-S test ta- 
ble is too conservative, i.e. with an appropriately generated 
look-up table (via Monte Carlo simulations), one can reject 
the null hypothesis that a sampled distribution is Normal at 
a higher significance level than with the standard table. 

The implications of this to our case are not known, 
and have not been investigated. Speculatively extrapolat- 
ing Lilliefors's conclusion, the 2D K-S test for comparing 
an unknown to a known distribution should be, if anything, 
more powerful than presented. This would increase the sig- 
nificance of source X6, making its association with the sug- 
gested BY Dra variable more likely. 

The advantage of an ML approach here would be that 
likelihood ratios between different models do not suffer from 
such problems. 



7.2 Significance of the detections 

The developed detection significance test for additional 
may seem subjective, since prior 



sources in Section 4.2 



knowledge is needed about the Z„ • curve (where M is the 
number of sources in the image). Naturally, this information 
is not available when working with an astronomical image 
representing an unknown source distribution, where M is a 
sought parameter. However, due the (nearly) distribution- 
free character of the pseudo one-sample 2D K-S test (Sec- 
tion 3.1.2), all that is needed is the correlation coefficient 
CC of the counts in the image, which is readily available 
(CC = 0.17 for the ROSAT image in Figure |). Provided 
that the best-fit model with TV sources represents the im- 
age reasonably well (K-S probability > 1 per cent), the 
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curve will be indistinguishable from the Z„ ' curve 
of the image, since the correlation coefficients of the A'^- 
source model and of the (M-source) image will be very sim- 
ilar. Indeed, in our case the best-fit five-source model has 
CC = 0.19, which given the slow dependence of the Zn dis- 
tribution on CC, well approximates the Zn distribution for 
CC = 0.17 (the correlation coefficient of the counts in the 
ROSAT image). 

More general than the false-detection probability is the 
fraction of cases, in which the observed image can be 
represented by a best-fitting model containing K < N ^ M 
sources, where K does not necessarily equal M — 1. Here A'^ 
is our best guess for the number of sources in the image, and 
M is the actual (unknown) number of sources. The quan- 
tity 1 — P'jf is the significance level at which we can reject 



the hypothesis that the image contains only K sources. To 
determine using the 2D K-S test, we need to compare 
multiple images of the same field to a single model simula- 
tion with K sources (Section 3.1.2). Naturally, this cannot be 
done, since there rarely exist multiple available images of the 
same field. However, is well approximated by the quan- 
tity Pk,n, given that, as discussed above, Z^'^ describes 
well the Z^'''^' distribution of the image. This is the value 
listed in Table ^ (using N = 5) for the probability that the 
ROSAT image can be fitted with fewer than 5 sources. For 
K = N — 1, in Pk,n we recover the false-detection probabil- 
ity for the A''th source, as already discussed in Section 



4.2 



7.3 Detected sources 

The positions and the count rates of the detected sources are 
in excelle nt agreement with Model I (b ased on 1995 ROSAT 
data) of Verbunt fc Johnston (200C ) in their maximum- 



likelihood analysis of the same H RI field. A Chandra im age 
of the core region of NGC 6793 ( [Grindlay et al. 2001a| ) re- 
veals a greater complexity of sources (Figure ^) . The source 
"doubles" CVl and CV4, as well as CV2 and U18 are too 
close (~ 2.5 arcsec « 0.3 FWHM) to be distinguished as 
separate sources in the ROSAT/'RYU image, and are rep- 
resented as blended sources XI and X4, respectively. The 
remainder of the sources marked with open squa res are too 



faint (L^ < lO"*^ erg s"^ Grindlay et al. 2001a) to be de 



tected given the crowdedness of the field. None the less, there 
is a clear one-to-one correspondence between the brightest 
(L^ > 10^^ erg s"^) Chandra sources, and the ones detected 
in the ROSAT/YiRl image usin g the 2D K-S techniq ue. 



Although in their Model IV Verbunt & Johnston predict 
the existence of separate x-ray counterparts to sources CV2 
and U18, that model is fit to 1991 ROSAT/YiRl data when 
CV2 was more prominent in x-rays relative to U 18 (hence 
could be more accurately centroided; cf. Fig. 1 in Cool et 



al. 1993), and three of the sources in the model have fixed 
positions. On the other hand, in our 2D K-S test iterative 
analysis we have not used any fixed parameters. Moreover, 
the K-S test suggests the existence of source U43, detected 
(albeit inconclusively, and offset by ~ 4.5 arcsec from its 
Chandra position) as source X6, for which there exists no 
x-ray identifi cation prior to the Chandra results of Grindlay 
et al. (2001c). Althou gh the source is fainter (logL^, = 29.4 



Grindlay et al. 20Qla) than other undetected sources in the 
complex X1-X4, the source must have been > lOx brighter 
to have been detectable with ROSAT (indeed, BY Dra bina- 
ries were disc overed in globular clus ters as faint and flaring 
x-ray sources, Grindlay et al. 2001b). 



8 CONCLUSION 



We have developed an application of the 2D K-S test (Pea- 



cock 1982 



Fasano ( 



: Franceschini 1987) in a source-detection 
algorithm for astronomical images. By employing the "sub- 
pixelization" technique on 3D astronomical images, we show 
that the 2D K-S test has greater power than the 3D K-S test 
- the intuitive choice for such images. We use Monte Carlo 
integration to determine the cumulative values of the pro- 
posed model distribution in all four quadrants around each 
count and, recognising the deviations that this incurs from 
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the derived Zn distributions, we provide our own reference 
tables for estimating the K-S probability. 

We devise an iterative source-modelhng routine that 
employs the K-S probability as a goodness of fit estimator, 
and can be used to find the optimum number, positions, and 
intensities of blended sources. We then apply the iteration 
scheme to a deep (75 ksec) ROSAT/YiRI exposure of the core 
region of NGC 6397 and find five blended sources, as well as 
a possible sixth one. The locations of the five brightest (and 
possible sixth) x-ray sources match closely (within the posi- 
tional error bars) the locations of probable CVs and BY Dra 



system s, disco vered with HST ( Cool et al. 1993 ^ 1995 ; Tay- 
lor et at 2001), and confirmed with Chandra (Grindlay et al 



2001a|!' The sixth source, X6, is a marginal detection with 
the 2D K-S technique and is l ikely identified with th e much 
fainter Chandra source, U43 (Grindlay et al. 2001a), which 



is in turn identified with a BY Dra binary, PC-4 ( Taylor et 
al. 200]|). 

Comparisons to other source-detection schemes 
(sliding-cell, wavelet detect, daophot/allstar, L-R de- 
convolution and ML techniques) applied to the same image 
demonstrate the superior power of our method in heavily 
crowded fields with low signal-to-noise. The example with 
the ROSAT/Yi'Rl deep field indicates that the proposed 
iterative source-modelling scheme can find applications in 
small-number statistics high-energy imaging, e.g. in deep 
exposures of globular clusters and extragalactic nuclear 
regions with Chandra, where the size of the PSF is often 
comparable to the angular separation between the objects. 

We thank our collaborators Adrienne Cool and Peter 
Edmonds for numerous discussions. This research was par- 
tially supported by NASA/LTSA grant NAG5-3256 and by 
HST grant GO-06742. 
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